function out = unconMoments(ySimAsInData,setupFilter)
FontSizeValue   = 14;
figureSize      =  [0 0 1000 600]; % (0,0) size in x size in y
numLags         = 20;
numBoots        = 10000;
bootWidth       = 60;     %all auto correlations in data are about zero at 50 lags

%% Computing moments in data and the model
nyMom         = setupFilter.nyMom;
T             = setupFilter.T;
numSim        = size(ySimAsInData,2);
numRep        = floor(numSim/T);

data          = setupFilter.data;
meanData      = zeros(nyMom,1);
stdData       = zeros(nyMom,1);
slopeData     = nan(nyMom,1);
crossData     = corr(data(:,1:nyMom));
autoCorrData  = zeros(nyMom,numLags);
skewData      = zeros(nyMom,1);
kurtData      = zeros(nyMom,1);

meanModel      = zeros(nyMom,1);
stdModel       = zeros(nyMom,1);
slopeModel     = nan(nyMom,1);
crossModel     = corr(ySimAsInData(1:nyMom,:)');
autoCorrModel  = zeros(nyMom,numLags);
autoCorrModel_tmp = zeros(nyMom,numLags,numRep);
skewModel      = zeros(nyMom,1);
kurtModel      = zeros(nyMom,1);

for i=1:nyMom
    % Data moments
    meanData(i,:)      = mean(data(:,i));
    stdData(i,:)       = std(data(:,i));
    if i > 6
        slopeData(i,1) = mean(data(:,i)-data(:,6)); %long-yield - 1-year yield
    end
    tmp                = autocorr(data(:,i),numLags);
    autoCorrData(i,:)  = tmp(2:end)';
    skewData(i,1)      = skewness(data(:,i));
    kurtData(i,1)      = kurtosis(data(:,i));
    
    % Model moments    
    meanModel(i,1)      = mean(ySimAsInData(i,:));
    stdModel(i,1)       = std(ySimAsInData(i,:));
    if i > 6
        slopeModel(i,1) = mean(ySimAsInData(i,:)-ySimAsInData(6,:)); %long-yield - 1-year yield
    end
    tmp                = autocorr(ySimAsInData(i,:),numLags);
    autoCorrModel(i,:) = tmp(2:end)';
    for s=1:numRep
        tmp            = autocorr(ySimAsInData(i,1+(s-1)*T:s*T),numLags);
        autoCorrModel_tmp(i,:,s) = tmp(2:end)';
    end
    skewModel(i,1)     = skewness(ySimAsInData(i,:));
    kurtModel(i,1)     = kurtosis(ySimAsInData(i,:));
end

%% Bootstrapping empirical moments
bsdata           = block_bootstrap(data(:,1:nyMom),bootWidth,numBoots);
autoCorrDataBoot = zeros(nyMom,numLags,numBoots);
meanDataBoot     = zeros(nyMom,numBoots);
stdDataBoot      = zeros(nyMom,numBoots);
slopeDataBoot    = nan(nyMom,numBoots);
skewDataBoot     = zeros(nyMom,numBoots);
kurtDataBoot     = zeros(nyMom,numBoots);
crossDataBoot    = zeros(nyMom,nyMom,numBoots);
for s=1:numBoots
    for i=1:nyMom
        tmp = autocorr(bsdata(:,i,s),numLags);
        autoCorrDataBoot(i,:,s) = tmp(2:end)';
        if i > 6
            slopeDataBoot(7:end,s) = mean(bsdata(:,7:end,s)-repmat(bsdata(:,6,s),1,4),1);
        end
    end
    meanDataBoot(:,s)    = mean(bsdata(:,:,s),1)';
    stdDataBoot(:,s)     = std(bsdata(:,:,s),0,1)';
    skewDataBoot(:,s)    = skewness(bsdata(:,:,s),0,1)';
    kurtDataBoot(:,s)    = kurtosis(bsdata(:,:,s),0,1)';
    crossDataBoot(:,:,s) = corr(bsdata(:,:,s));
end
meanDataCI95_prcMethod = prctile(meanDataBoot,[2.5 97.5],2);
meanDataSE        = std(meanDataBoot,0,2);
meanData          = mean(data(:,1:nyMom),1)';
meanDataCI95      = [meanData-1.96*meanDataSE,meanData+1.96*meanDataSE];

slopeDataCI95_prcMethod = prctile(slopeDataBoot,[2.5 97.5],2);
slopeDataSE      = std(slopeDataBoot,0,2);
slopeDataCI95    = [slopeData-1.96*slopeDataSE,slopeData+1.96*slopeDataSE];

stdDataCI95_prcMethod = prctile(stdDataBoot,[2.5 97.5],2);
stdDataSE        = std(stdDataBoot,0,2);
stdData          = std(data(:,1:nyMom),0,1)';
stdDataCI95      = [stdData-1.96*stdDataSE,stdData+1.96*stdDataSE];

autoCorrDataCI95_prcMethod = prctile(autoCorrDataBoot,[2.5 97.5],3);
autoCorrDataSE             = std(autoCorrDataBoot,0,3);
autoCorrDataCI95           = zeros(nyMom,numLags,2);
autoCorrDataCI95(:,:,1)    = autoCorrData-1.96*autoCorrDataSE;
autoCorrDataCI95(:,:,2)    = autoCorrData+1.96*autoCorrDataSE;

crossDataCI95_prcMethod = prctile(crossDataBoot,[2.5 97.5],3);
crossDataSE             = std(crossDataBoot,0,3);
crossDataCI95           = zeros(nyMom,nyMom,2);
crossDataCI95(:,:,1)    = crossData-1.96*crossDataSE;
crossDataCI95(:,:,2)    = crossData+1.96*crossDataSE;

skewDataCI95_prcMethod = prctile(skewDataBoot,[2.5 97.5],2);
skewDataSE       = std(skewDataBoot,0,2);
skewDataCI95     = [skewData-1.96*skewDataSE,skewData+1.96*skewDataSE];

kurtDataCI95_prcMethod = prctile(kurtDataBoot,[2.5 97.5],2);
kurtDataSE       = std(kurtDataBoot,0,2);
kurtDataCI95     = [kurtData-1.96*kurtDataSE,kurtData+1.96*kurtDataSE];


%% Plotting
figure('Name','Persistence','NumberTitle','off','Position',figureSize)
for i=1:nyMom
   subplot(4,3,i)
   hold on 
   tmp0 = plot(autoCorrData(i,:),'-k');
   tmp1 = plot(autoCorrData(i,:)-1.96*autoCorrDataSE(i,:),'--k');
   tmp2 = plot(autoCorrData(i,:)+1.96*autoCorrDataSE(i,:),'--k');
   tmp3 = plot(autoCorrModel(i,:),'-ok')  ;
   axis([1 numLags -.20 1.15])
   title(setupFilter.dataLabels{i},'Interpreter','Latex','FontSize',FontSizeValue+2)
   set(gca,'FontSize',FontSizeValue)
   if i == 1
       legend([tmp1(1), tmp0(1), tmp3(1)],{'CI 95 pct','Data','Model'},'Orientation','horizontal','FontSize',FontSizeValue+2,'EdgeColor',[0.15 0.15 0.15],...
           'Position',[0.384822896027229 0.939478465696058 0.243629479375424 0.0675996218505801],'FontSize',FontSizeValue);
       legend boxoff  
   end
end

%% Output
out.meanModel               = meanModel;
out.meanData                = meanData;
out.meanDataSE              = meanDataSE;
out.meanDataCI95            = meanDataCI95;
out.meanDataCI95_pcrMethod  = meanDataCI95_prcMethod;

out.slopeModel              = slopeModel;
out.slopeData               = slopeData;
out.slopeDataSE             = slopeDataSE;
out.slopeDataCI95           = slopeDataCI95;
out.slopeDataCI95_prcMethod = slopeDataCI95_prcMethod;

out.stdModel                = stdModel;
out.stdData                 = stdData;
out.stdDataSE               = stdDataSE;
out.stdDataCI95             = stdDataCI95;
out.stdDataCI95_prcMethod   = stdDataCI95_prcMethod;

out.autoCorrModel           = autoCorrModel;
out.autoCorrModelFinite     = mean(autoCorrModel_tmp,3);
out.autoCorrData            = autoCorrData;
out.autoCorrDataSE          = autoCorrDataSE;
out.autoCorrDataCI95        = autoCorrDataCI95;
out.autoCorrDataCI95_prcMethod = autoCorrDataCI95_prcMethod;

out.crossModel              = crossModel;
out.crossData               = crossData;
out.crossDataSE             = crossDataSE;
out.crossDataCI95           = crossDataCI95;
out.crossDataCI95_prcMethod = crossDataCI95_prcMethod;

out.skewModel               = skewModel;
out.skewData                = skewData;
out.skewDataSE              = skewDataSE;
out.skewDataCI95            = skewDataCI95;
out.skewDataCI95_prcMethod  = skewDataCI95_prcMethod;

out.kurtModel               = kurtModel;
out.kurtData                = kurtData;
out.kurtDataSE              = kurtDataSE;
out.kurtDataCI95            = kurtDataCI95;
out.kurtDataCI95_prcMethod  = kurtDataCI95_prcMethod;

out.numBoots         = numBoots;
out.bootWidth        = bootWidth;


end